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Abstract. Iteratively Re-weighted Least Squares (IRLS) were used to simulate the 
Lp-norm approximation of the ballistic trajectory in absolute gravimeters. Two 
iterations of the IRLS delivered sufficient accuracy of the approximation without 
a signihcant bias. The simulations were performed on different samplings and 
perturbations of the trajectory. For the platykurtic distributions of the perturbations, 
the Lp-approximation with 3 < p < 4 was found to yield several times more precise 
gravity estimates compared to the standard least-squares. The simulation results were 
confirmed by processing real gravity observations performed at the excessive noise 
conditions. 
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1. Introduction 

Absolute ballistic gravimeters measure gravity acceleration by tracking the free motion 
of the test mass in the gravity held. Having the distances {Si,S n} covered by the 
test mass over the time intervals {Ti,..., T/v}, the acceleration is found as parameter 
of some trajectory model htted to the data pairs (Tj, Si) by methods of the regression 
analysis. One of the commonly used models describes unperturbed vertical motion of 
the test mass in the gravity held 

Zi = zq + VoTi + g /2. 


( 1 ) 
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The expression (1) represents a linear regression model, because it’s defined by a linear 
combination of the estimated parameters Zq, The estimates are usually found by 

minimising the sum of the squared discrepancies between the measured data and the 
model: 

{zo, V^o, g) : X] e- ^ min, where ei = Si - Zi. (2) 

Minimised here is the square norm (L 2 -norm) of the vector e = {e*}, which is a special 
case of the more general Lp-norm 

l|eilp = (5^ l^r) P > 1- (3) 

The tradition of minimising the sum of squares (i.e. finding the least squares) has a long 
and reach history dating back to the beginning of the 19-th century. At that time the 
main reason for using the least squares was their computational simplicity, as for the 
p=2 the fitting of a linear model is reduced to finding a unique solution for a full-rank 
system of linear equations. At any other p the equations become non-linear in estimated 
parameters and require involved computations to find the solution. 

Another reason for using the least squares relates to a different kind of linearity: 
linearity of the estimates with respect to the dependent variable {Si in our case), which 
implies no bias of the estimates. Moreover, according to the Gauss-Markov theorem, of 
all the linear/unbiased estimates, the least-squares deliver the most precise one, in case 
the errors in the observations are uncorrelated, have zero mean and constant variance. 
This is especially important in metrology, where the bias may lead to an incorrect 
estimates of the measurand and its uncertainty. 

Nowadays, when advanced computing is no longer a problem, the Lp-approximation 
with p 7 ^ 2 is widely used [l]-[7]. Many researches have noticed that the value of p 
delivering the lowest variance to the Lp-estimates decreases as the kurtosisf of the error 
distribution grows; there are formulae suggested for this dependence [8, 9]. Most results 
have been obtained for the values 1 < p < 2, that provide better estimates when the 
noises have heavily-tailed distributions. In trajectory tracking, however, more prevalent 
are limited-band noises, for which the properties of the Lp-approximation are less known. 

The purpose of this work is to investigate the feasibility and possible advantages of 
approximating the gravimeter’s test mass trajectory in different Lp-norms. 

The paper has the following structure. The chapter 2 discusses implementation 
of the Iteratively Re-weighted Least Squares (IRLS) for the Lp-norm approximation. 
The chapter 3 describes the Monte-Carlo simulation of gravity estimates in Lp-norm. A 
real-life example of improving the estimates by the Lp-norm approximation is discussed 
in the chapter 4. We summarise our study in the chapter 5. 

f The kurtosis |32 characterises the spread of the distribution over the possible values of the random 
variable. It’s defined as pn /where p .4 is the fourth central moment, cr is the standard deviation. 
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2. The Lp-norm approximation by the Iteratively Re-weighted least-sqnares 


We implement the Lp-norm approximation by performing the weighted least-squares 
(WLS) approximation, in which the weights are dehned in a special way [10]. For the 
model (1), the WLS solution is determined by 

« = A+S, (4) 

where 

^ = {zq, Vo,g)'^ is the vector of the unknown parameters, 

S = {Si, 5*2,..., Lat)^ is the vector of the observations, 

A+ is the Moore-Penrose pseudo-inverse matrix calculated as 

A+ = {A^VA)-^A^V, (5) 

where 
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Till \ 
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n/2) 



is the experiment design matrix. 


W = diag(ui,U 2 , •••Wv) (7) 

is the diagonal matrix of weights. The approximation in the Lp-norm is achieved when 
the weights Vi are related to the discrepancies e* like [10] 

Vi = ( 8 ) 

For p=2, the weights (8) turn into units, yielding the standard least-squares solution. 
The formula (8) cannot be used directly to find the weights Vi, because the weights 
need to be known before the approximation, while the residuals e* are known only after 
the approximation is done. Because of that, the WLS is usually applied several times 
leading to the Iteratively Re-weighted Least Squares. The iterations start with the 
weights Vi=l, i.e. with the regular least-squares. The residuals e* are then used to build 
the new weights 

V = diag(|ei|P“^ |e2^■^ ..., (9) 


which are then applied to approximate the original data. The iterations can be repeated, 
but the process is known to sometimes have a problematic convergence. Though 
different approaches exist to improve the convergence [10], the simulations found that 
two iterations always produce a sufficiently accurate for our purposes approximation, 
with the bias not discernible in random noise. For better computational stability we 
normalised the weights (8), so they would range from 0 to 1: 


Vi = 


M ) 


L > 2, 


( 10 ) 
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M = max I ei 


( 11 ) 


If p < 2, the powers in (8) become negative, and to prevent zero division errors another 
modihcation of the weights was necessary: 

max( I ei I, rM)y“^ 

Vm ) ’ ^ ^ ’ 




where r = 0.001 is a regularisation parameter establishing the threshold of smallness 
as fraction r of the maximum residual M. According to (12), the residuals below the 
threshold get replaced with the threshold value rM. A similar way of regularisation is 
implemented in [7]. 


3. Modelling the properties of the absolute gravity estimates in the 
Lp-norm 

3.1. Organisation of the simulation experiment 

For the simulations, we added random perturbations to the parabolic trajectory (1) 
with known acceleration g = 9.8 ms“^§, and then recovered the acceleration from the 
perturbed parabola using the Lp-norm approximation (hgure 1). The values of p ran 
from 1 to 6 with the increment of 0.1, so the standard least-squares solution (p=2) was 
also available from the simulation and was used for comparison. To build the parabola, 
700 data points {Ti, ...,T 7 oo} representing 699 time intervals with the common start at 
Ti were used, spanning the total time T 700 = 0.22 s. We considered 3 distributions 
of the Tj-points (experiment designs) found in different types of absolute gravimeters 
(hgure 2). The EST (equally spaced in time) design has the uniform distribution of the 
Tj-points and is used in both free-fall and rise-and-fall gravimeters. The ESD^p design 
has the minimal density at the start, linearly increasing towards the end of the interval. 
The ESDj^p design has the minimal density in the centre, linearly increasing towards the 
start and the end of the interval. These two designs are found in free-fall (FF) and rise- 
and-fall (RF) gravimeters with the levels equally spaced in distance (ESD). For every 
type of perturbation we simulated 3000 drops and, after approximating the trajectory 
in different Lp-norms, analysed the results. Every realisation of the simulated noise 
was used in all three experiment designs by assigning either EST, ESDpp, or ESD^^^p 
realisations of the vector T = {Ti, ...,T 7 oo} to the same noise vector e = {ei, ..., € 700 }- 
The simulations and data analysis were performed with the MATLAB® software. 

3.2. Simulation results for the random uncorrelated perturbations 

The trajectory noise in absolute gravimeters comes from a number of sources, such as 
ground vibrations, laser stabilisation, digital signal acquisition, fringe counting, etc., so 

§ Other parameters were Zq = 1 nm and Vq = 0 ms“^. Zero initial velocity was chosen to maximise 
the difference between the experiment designs. 
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Figure 1. Numerical simulation of the properties of the ballistic trajectory 
approximation in the Lp-norm. 





Figure 2. Experiment designs used in the simulation: 

EST - equally spaced in time. Time intervals do not depend on the trajectory, can be 
used in both free-fall and rise-and-fall gravimeters; 

ESDpp - equally spaced in distance for the free-fall trajectory; 

ESDpp - equally spaced in distance for the rise-and-fall trajectory. 


the combined noise most often is of the Gaussian type. Sometimes, however, individual 
sources can dominate, creating other noise distributions. We did the simulations for 
several noise types, all having the standard deviation of 1 nm. 
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For all noise distributions, the distribution of the gf-estimates was normal, for any 
p. For every noise distribution there existed an optimal value of p delivering minimal 





P 


Figure 3. Three ways of graphical representation of the simulation results (noise: 
Uniform, trajectory sampling: EST) 

a - mean gravity estimates shown relatively to the known value oi g = 9.8 ms“^ (1 
pGal = 10“® ms“^) and their standard deviations as error bars, for different L^-norms; 
b - variance of the Lp-norm estimates relative to the variance of the standard least- 
squares solution; 

c - relative efficiency of the estimates as ratio of the variance of the L 2 -estimate to the 
variance of the Lp-estimate - shows the number of times the variance has decreased 
compared to the standard least-squares. 


variance of the estimate (hgure 3a). The optimal p ranged from 1.4 to 3.3 depending on 
the noise distribution. This result deviates from the theory [11] predicting no optimal 
p for at least the Uniform distribution. 

For the numeric characteristic of the distribution shape, we used anitkurtosis||, 
because it can only assume values from 0 to 1 for any distribution. We found it more 

II The antikurtosis x relates to the kurtosis (32 as x = l/vlU [12]. 
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convenient than the traditionally used kurtosis, which can range from 1 to infinity. 
The optimal p increased with antikurtosis (table 1), but not always the variance at 
the optimal p improved substantially. We visualised the improvement by plotting the 
variance for every p (figure 3b) divided by the variance of the standard least-squares. The 
factor of the variance improvement is shown on the plot of the inverse ratio (figure 3c), 
which is equivalent to the relative efficiency of the Lp-estimate compared to the standard 
least-squares one. The relative efficiency plots for all simulated noises and experiment 
designs are presented on the figure 4. As seen on the figure and at the table 1 , significant 
improvement of the precision is achieved for the Uniform and Arcsine distributions, if 
the norm L 3.3 is used for the parabola fitting instead of the L 2 . 

As non-linearity of the Lp-estimates can entail bias, we did the drop-to-drop 
comparison of the L 2 -estimates (theoretically unbiased) with the Lp-estimates at the 
optimal p. No statistically significant divergence of the two estimates was detected. 
We thus confirmed no bias for our approximation procedure, the fact known for the 
Lp-estimates in case of the symmetric noise distributions [ 8 , 11]. 



Figure 4. Relative efficiency of the Lp-norm estimates (Y-axis) for different noise 
types. The efficiency is calculated as ratio of the L 2 -estimate variance for the EST- 
design to other Lp-estimates for all designs. 
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Table 1. Results of the simulation of the absolute gravity estimates by approximating 
the trajectory in different Lp- norms. The case of random uncorrelated perturbations 
of the trajectory. 


Noise 

Distribution 

Anti¬ 

kurtosis 

X 

Best 

P 

Rel. 

eff.®' for the best p 

EST 

ESDpp 

ESD^p 

Laplace 

0.41 

1.4 

1.3 

1.3 

1.3 

Normal 

0.58 

2.0 

1.0 

1.0 

1.0 

Triangle 

0.65 

2.5 

1.1 

1.1 

1.1 

Uniform 

0.75 

3.3 

2.8 

2.8 

2.8 

Arcsin 

0.82 

3.3 

7.0 

7.0 

7.1 


^'Determined as ratio of the variance of the L 2 -estimate to the variance of the Lp- 
estimate for the same experiment design. 


3.3. Simulation results for the “harmonie in the Gaussian noise” perturbations 

According to the table 1, the Lp-estimates with p > 2 can provide a signihcant gain 
in the efficiency when the errors are distributed by the Uniform or Arcsine law. In 
trajectory tracking such errors are often caused by harmonic perturbations of the 
reference frame. We modelled the perturbations by the sum of a sinusoid and Gaussian 
noise. In every simulated drop, a sinusoid with a given frequency / and amplitude of 
1.41 nm was assigned an initial phase randomly taken from the uniform distribution 
on [0,27r], and then sampled according to the experiment design. The sinusoid samples 
were added by the Gaussian noise, one noise realisation for all three experiment designs. 
The signal-to-noise ratio (SNR) is dehned as quotient of the standard deviations of the 
sinusoid and the noise, so the inhnite SNR means no random noise. 

The results of the simulation for the frequencies of 17 Hz, 35 Hz, 55 Hz and the SNR 
values of oo, 100, 10, 2 are shown in the table 2. The antikurtosis x of the perturbation 
was from 0.70 to 0.82, typical for the range occupied by the Uniform and Arcsine 
distributions. The resulting gf-estimates have passed the normality test (Lilliefors , 5%) 
only for some low SNR values. The optimal p and relative efficiencies for which the 
normality was observed are shaded in the table 2. Like for uncorrelated noises (table 1), 
there always existed an optimal p delivering the lowest value to the variance of the 
Lp-estimate. The optimal p was between 2.8 and 3.8, yielding the variance decrease of 
the estimate between 1.2 and 30.6 times. The values of p and the efficiency gain have a 
complicated dependence on the signal frequency, the experiment design, and the SNR, 
so the values have to be evaluated experimentally in every practical case. 

4. Example of using the Lp-estimates on real data 

The held measurement of absolute gravity are often performed under the harsh 
geophysical conditions that deteriorate the performance of the instruments. Oftentimes 
it’s not possible to repeat the measurements or improve the conditions, so the 
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Table 2. Results of the simulation of the absolute gravity estimates by approximating 
the trajectory in different Lp- norms. The trajectory perturbations are “harmonic in 
the Gaussian noise.” 


//Hz 

SNR 


EST 


ESDpp 


ESD^, 

X 

P 

rel.eff.®' 

X 

P 

rel.eff.®' 

X 

P 

rel.eff.®' 

17 

oo 

0.81 

3.6 

5.5 

0.81 

3.5 

17.5 

0.78 

3.4 

1.5 


100 

0.81 

3.6 

5.8 

0.81 

3.5 

16.9 

0.78 

3.4 

1.5 


10 

0.80 

3.6 

5.5 

0.80 

3.5 

17.7 

0.78 

3.4 

1.6 


2 

0.70 

3.6 

2.7 

0.70 

3.8 

7.4 

0.68 

3.3 

1.2 

35 

oo 

0.82 

3.5 

29.3 

0.82 

3.4 

22.2 

0.82 

3.6 

14.9 


100 

0.82 

3.5 

29.5 

0.82 

3.4 

23.2 

0.82 

3.5 

15.1 


10 

0.81 

3.5 

22.0 

0.81 

3.4 

19.6 

0.81 

3.6 

13.7 


2 

0.70 

3.5 

2.7 

0.70 

3.4 

3.3 

0.70 

3.5 

3.3 

55 

oo 

0.82 

3.3 

16.0 

0.82 

3.3 

28.5 

0.82 

3.3 

25.8 


100 

0.82 

3.3 

15.8 

0.82 

3.3 

28.8 

0.82 

3.3 

25.6 


10 

0.80 

3.3 

15.6 

0.80 

3.3 

30.6 

0.81 

3.3 

25.3 


2 

0.70 

2.8 

1.3 

0.70 

3.1 

2.3 

0.70 

3.1 

2.2 


^'Determined as ratio of the variance of the L2-estimate to the variance of the Lp- 
estimate for the same experiment design. 


researches have to make the best possible gravity estimates using the available data. 
As example, we consider a series of measurements taken by the TAG-1 absolute 
gravimeter developed by the Earthquake Research Institute of the University of Tokyo. 
The instrument features a laser interferometer with built-in accelerometer for the 
independent measurement and correction of the ground vibrations, and a homodyne 
quadrature fringe signal detector [13, 14]. The free-fall trajectory of about 12 cm is 
sampled with the frequency of 10 MHz (EST-design). A typical drop-to-drop scatter of 
the gf-values at the station ‘Esashi’ is about 20 pGal^, so the standard deviation of the 
mean below 1 pGal is reached in several hundreds drops. The bandwidth of the built-in 
accelerometer is up to 4 Hz, so the increased intensity of the high-frequency seismic noise 
could signihcantly disturb the measurements. The hgure 5(a) shows one such noisy set 
of gravity estimates taken on December 23, 2011 at the ‘Esashi’ gravity station, after 
applying standard tidal, atmospheric, and polar motion corrections. To obtain the 
Lp-estimates, the trajectory was decimated to 2 kHz, the standard speed-of-light and 
gradient corrections were applied to the time and distance coordinates respectively, 
then the gravity was found using the standard least-squares approximation. The set 
fails the normality test, making the standard ‘three-sigma’ rejection too restrictive for 
the obtained heavily-tailed distribution of the L 2 -estimates. However, neither four, nor 
hve sigmas made much difference to the mean g-value or its standard deviation, which 
remained at about 3 p.Gal (table 3). 

We then approximated the trajectory in the Lp-norm for p ranging from 1 to 6 (see 


^ 1 pGal =10 ® ms ^ 
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Figure 5. Improving precision of the absolute gravity estimates by using the Lp-norm 
approximation of the free-fall trajectory: 

(a) standard least-squares approximation (p=2), the bars are for the 3ff, 4cr, and 5cr; 

(b) approximation for p=3.5, the bars are for the 3cr, 4cr, and 5cr; 

(c) mixed p=2/p=3.5 approximation based on the antikurtosis of the residuals: black 
dots - Gaussian group [p = 2), grey dots - Arcsine group (p = 3.5). 


section 2). For each p we obtained the variance of the set gravity estimate and plotted 
the relevant relative efficiency (hgure 6). According to the plot (hgure 6(a), curve 1), the 
variance at p=3.5 decreased several times w.r.t. p=2, reducing the standard deviation 
of the mean to about 1 pGal (table 3). Based on the simulations (section 3), the peak 
of the efficiency at p > 2 is suggestive of the residuals with non-Gaussian distributions. 
We investigated the residuals and found the values of antikurtosis ranging from 0.487 to 
0.813, which corresponds to the wide range of distributions - from Laplace to Arcsine 
(hgure 7). The likely reason for the non-Gaussian errors was resonant oscillations of the 
reference mirror caused by the excessive seismic noise. The most prevalent in the set 
were residuals with the Gaussian and Arcsine distributions. According to the simulations 
(section 3), the most efficient Lp-approximations for the Arcsine-noised trajectory were 
achieved for the values of p between 3 and 4, which explains the decrease of the 
uncertainty. However, for the Gaussian noises the approximation with p in that range 
can be several times less efficient than the standard least-squares (hgure 4, second from 
the top). This led us to the idea of splitting the set into two parts with predominantly 
Gaussian and Arcsine noises, and independent optimising the Lp-norm on each set. We 
have split the groups at the antikurtosis x = 0.65, about in the middle of two peaks at 
the hgure 7. Gravity estimates in each group have passed the normality test (Lilliefors 
, 5%), so further processing was carried out with the 3 ct rejection, each group with its 
own ff. The splitting has increased the efficiency of the estimate for the Arcsine group 
with p=3.5 from 10.0 to 12.5 (hgure 6(a), 2). For the Gaussian group we used the 
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p p 

Figure 6. Relative efficiency of the gravity estimates with the Lp-norm approximation 
of the trajectory: 

(a), 1: The entire set with the 5cr censoring. The variance of the Lp-approximation for 
p=3.5 decreases 10 times. 

(a) , 2: The portion of the set with high antikurtosis values. The variance of the Lp- 
approximation for p=3.5 decreases 12.5 times. 3cr censoring applied. 

(b) : The Gaussian portion of the set. The standard least-squares (p = 2) are close to 
optimal. 3cr censoring applied. 


Table 3. Estimates of the absolute gravity value obtained with different modes of the 
approximation of the free-fall trajectory 


Mode of approximation 
of the trajectory 

Rejection 

limits 

Accepted 

drops 

Set mean 

9 / M-Gal 

Set std.dev. Mean’s std.dev. 

(f / p.Gal ff / p.Gal 

Whole set, L 2 

3(T 

784 

980 121 363.52 

80.44 

2.87 


4cr 

803 

980 121 362.38 

95.24 

3.36 


5cr 

811 

980 121 363.96 

107.31 

3.77 

Whole set, L 3 5 

3(T 

802 

980 121 362.70 

29.51 

1.04 


4cr 

817 

980 121 363.31 

32.67 

1.14 


5cr 

820 

980 121 362.76 

33.88 

1.18 

Split groups, 12 / 13,5 

3o-/3ff 

396/406 

980 121 364.50’" 

16.37/38.18 

0.76'" 


^ Weighted average of two groups, the weights are the reciprocals of the groups’ vari¬ 
ances 

Standard deviation of the weighted average 


standard least-squares with p= 2 . 0 , which helped to avoid the decrease of the efficiency 
at p=3.5 for that group (hgure 5(b)). Combining the estimates for both groups has 
reduced the standard deviation of the weighted mean to 0.76 pGal (hgure 5(c)). The 
gravity estimates obtained with different approximations are summarised in the table 3. 
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Drops 


Figure 7. Antikurtosis composition of the noise in the data set. 

The histogram on the right-hand side shows the distribution of the antikurtosis for the 
noises found in the set’s drops. The histograms in the middle and the plots on the left- 
hand side show examples of the noises present in the set. Due to the inclement seismic 
conditions, the noises cover a wide range of distributions, from Laplace to Arcsine. 

5. Conclusions 

We have simulated the approximation of the ballistic trajectory in different Lp- norms, 
in view of using it to improve the absolute gravity estimates. To perform the 
approximation, the Iteratively Re-weighted Least Squares method was implemented. 
The estimates were simulated for three different plans of the experiment corresponding 
to two types of the ballistic trajectory (free-fall, rise-and-fall) and two schemes of data 
location (equally spaced in time or in distance). Each plan was applied to several 
perturbation types, including the cases of random uncorrelated noises and a harmonic 
in the Gaussian noise. As an application example, the simulation results were used to 
improve the absolute gravity estimates obtained at the ‘Esashi’ gravity station by the 
TAG-1 gravimeter. 

The following conclusions can be drawn from the study. 

(i) The Iteratively Re-weighted least-squares with two iterations approximate the 
ballistic trajectory in the Lp-norm with the accuracy sufficient for the purposes 
of absolute gravimetry. No statistically significant bias of the estimate was found 
for the wide range of the trajectory noises. 

(ii) The Lp-norm approximation of the ballistic trajectory with p > 2 estimates 
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gravity more precisely than the standard least-sqnares (p=2), if the trajectory 
pertnrbations deviate from the Gaussian towards more platykurtic distributions, 
i.e. those with flattish or concave shape. 

(iii) The observed efficiency maximum for certain values of p on platykurtic distributions 
disagrees with the prediction of the work [11], where the monotonous increase of 
the efficiency with the norm power p is obtained for the Uniform distribution. This 
could possibly be explained by a limited iterations of our IRLS procedure, resulting 
in a “quasi” rather than the “exact” Lp-approximation. 

(iv) For the same measurement interval T and hxed number of data points iV, the 
EST sampling of the free-fall trajectory makes the gravity estimates about twice as 
efficient as those with the data equally spaced in distance, for any p. This points 
to an obvious way of improving the modern free-fall absolute gravimeters, which 
predominantly use the ESD scheme. For the symmetric rise-and-fall trajectory, the 
efficiency of the EST and ESD sampling is about the same. 

The presented study is a hrst attempt to obtain the gravity estimates by approximating 

the ballistic trajectory in the Lp-norm with p ^ 2. Many important questions are yet 

to be answered by future work. The questions include, but are not limited to 

(i) The scope of the applicability of the Lp-norm approximation. For the purposes of 
absolute gravimetry we investigated only the estimates of the quadratic coefficient 
of the second-order polynomial model. The results can be applied to other models 
to build more efficient estimates in wider range of applications. 

(ii) Auto-correlation of the residuals. The main source of the platykurtic distributions 
is the low-frequency noises. Our simulations of the correlated Arcsine noises and 
processing of the real data suggest that the efficiency of the Lp-approximation 
is dehned mostly by the noise distribution and is not deteriorated by the auto¬ 
correlation of the residuals. Still, the influence of the auto-correlation calls for 
more thorough investigation, including the comparison to other algorithms of the 
low-frequency noise hltering. 

(iii) Contaminated distributions within a drop. Further investigations are desirable to 
better understand the properties of the Lp-norm approximation for the mixed noise 
distributions, such as harmonics in the Gaussian noise. 

(iv) Contaminated distributions of the set and set splitting. When a set has drops with 
different noise types, a single value of p may not yield the best Lp-approximations 
for all drops. In the example (section 4) we have split the set in two at the centre of 
the antikurtosis histogram. Though we gained certain uncertainty decrease, there 
may exist better ways of set splitting, worth to be investigated. 

(v) The reasons of the departure of the simulated results from those predicted in [11] 
have to be investigated to better understand the properties of the IRLS algorithm 
of the Lp-approximation. 
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